Pearson Correlation RNA-seq only

Introduction

This is global analysis of all of the muations for all patients that we have RNA-seq for.

Only variant with at least 3% AF and 100x coverage at the given location (must be in all samples from all patient) will be kept. I will also remove patients whose bulk RNA-seq clusters in a weird way (e.g. GC RNA-seq does not cluster with other GC samples).

For any position that is not fulfiling the above criteria the position is dropped.

Finally, I remove troublesome MT position: 295, 2617, 13710 - reported as mutated in RNA-seq 3107 - This position is next to an deletion and RNA-seq and WGS do not work well together. 302 - 317 - This is the variable region of d-loop. When I manually check it there are lot of mutations in a stretch (C)7T(C)5 and they affect alignment of T nucleotide in a stretch https://www.nature.com/articles/cr200969

Setup

Nucleotide frequency

The following code can be run in bash to collect the frequencies of nucleotides. Requires samtools (http://www.htslib.org/) and pysamstats (https://github.com/alimanfoo/pysamstats). Here I used samtools 1.9 and pysamstats 1.1.2 (pysam 0.15.2)


for i in 0062 0231 0412 0553 0832 0848 0931 1075 1218 1281 1533 1546 1566 1661 1790;

do
for ii in BE GC NE D2 ;
  do
    echo $i
    echo $ii
#   samtools index ../AHM_$i"_"$ii"_Aligned_MT.bam"
############ RNA-seq data only have 3 qualities or read aligment based on the number of regions the read was mapped to. Keep only reads that map to a single location (-min-mapq=200)
############ The quality of based was variable and I only keep the best nucleotides with base quality above 34
    pysamstats --type variation --min-mapq=200 --min-baseq=34 -f ../../Files/hg37/Homo_sapiens.GRCh37.dna.chromosome.MT.fa ../../BAM/RNA-seq/AHM_$i"_"$ii"_Aligned_MT.bam" >./AHM_$i"_"$ii"_Aligned_MT_q.counts"
  done
done

Functions

library("edgeR")
library("ggplot2")
library("pheatmap")
library("tidyr")
library("reshape2")
library("RColorBrewer")
library("corrplot")

# function for calculation of mitodistance between samples
mitodist <- function(x, y, z = 100, vaf = 0.03, count = 0, positions = NA) {
    
    # z - coverage vaf - minimum VAF count - combined count in two samples for given
    # VAF
    
    # merge locations as they do not always overlap
    tmp <- merge(x, y, by = c("chrom", "pos"), all = TRUE, sort = TRUE)
    tmp <- tmp[!(tmp$pos %in% positions), ]
    # replace the NA with 0 at the counts table
    tmp[, c(4, 7:10, 12, 15:18)][is.na(tmp[, c(4, 7:10, 12, 15:18)])] <- 0
    
    # make matrix of frequencies
    x2 <- as.matrix(tmp[, 7:10]/rowSums(as.matrix(tmp[, 7:10])))
    x2[is.na(x2)] <- 0
    y2 <- as.matrix(tmp[, 15:18]/rowSums(as.matrix(tmp[, 15:18])))
    y2[is.na(y2)] <- 0
    # change counts to 0 if they are below indicated vaf in both samples (1 sample
    # above given vaf is required to keep the allel)
    tmp[, 7:10][x2 < vaf & y2 < vaf] <- 0
    tmp[, 15:18][x2 < vaf & y2 < vaf] <- 0
    
    # change counts to 0 if they are below required coverage for the given alelle
    # (default is minimum of 5 reads per pair/sum of samples)
    tmp[, 7:10][(tmp[, 7:10] + tmp[, 15:18]) < count] <- 0
    tmp[, 15:18][(tmp[, 7:10] + tmp[, 15:18]) < count] <- 0
    
    # recalculate the frequencies
    x2 <- as.matrix(tmp[, 7:10]/rowSums(as.matrix(tmp[, 7:10])))
    x2[is.na(x2)] <- 0
    y2 <- as.matrix(tmp[, 15:18]/rowSums(as.matrix(tmp[, 15:18])))
    y2[is.na(y2)] <- 0
    
    # get the indicator tables. This table check the coverage at each position. It
    # has to be larger than indicated coverage (default is 100)
    x1 <- as.matrix(ifelse(rowSums(as.matrix(tmp[, 7:10])) > z, 1, 0))
    y1 <- as.matrix(ifelse(rowSums(as.matrix(tmp[, 15:18])) > z, 1, 0))
    
    
    # calculate square root of absolute AF difference
    AF_dist <- sqrt(abs(x2 - y2))
    
    # calculate distance
    dist <- (sum(AF_dist * as.vector(x1 * y1)))/sum(x1 * y1)
    
    return(dist)
}

Check phenotypes

First, I will check the phenotypes of the samples to confirm that they are what they should be. In another analysis I noticed that some samples seem to cluster with wrong tissue, i.e. NE clusters with gastric tissue.

All Samples

# Names of samples:
samples <- c("0062", "0231", "0412", "0553", "0832", "0848", "0931", "1075", "1218", 
    "1281", "1533", "1566", "1661", "1790")  #RNA-seq data
conditions <- c("BE", "GC", "D2", "NE")  # RNA-seq data

files.dir <- "~/Dropbox/Postdoc/2019-12-29_BE2020/RNA-seq/"

# EAC and normal data (tmp, linear scale)
bulk.data <- readRDS(paste0(files.dir, "/EAC-all/count_codingGenes.rds"))
IDs <- grep(paste(samples, collapse = "|"), colnames(bulk.data), value = TRUE)
bulk.data <- bulk.data[, IDs]

# get information about the samples from EAC
metadata <- readRDS(paste0(files.dir, "//EAC-all/samplesAnnotation.rds"))
metadata <- metadata[IDs, 6:8]
group <- factor(metadata$sampleGroup)  #Levels: Barretts Duodenum Gastric Squamous


y <- DGEList(bulk.data, group = group)

# #keep only genes with expression above 1 cpm in at least 3 samples
keep <- rowSums(cpm(y) > 1) >= 3
y <- y[keep, , keep.lib.sizes = FALSE]

# get normalization factors
y <- calcNormFactors(y)


# Start performing diff expression - get sample desing model design <-
# model.matrix(~ batch + group)
design <- model.matrix(~0 + group)
rownames(design) <- colnames(y)

colnames(design) <- levels(group)

# Estimate disparsion
y <- estimateDisp(y, design, robust = TRUE)

plotBCV(y)

# Fit general linear model data
fit <- glmQLFit(y, design, robust = TRUE)

plotQLDisp(fit)

points <- rep(c(0, 1, 2), 2)
colors <- rep(c("blue", "darkgreen", "red"), 2)


p <- plotMDS(y, col = colors[group], pch = points[group], plot = FALSE)


p2 <- data.frame(p$x, p$y, .Tissue = metadata$sampleGroup, .Batch = metadata$batch)


ggplot(data = p2) + geom_point(mapping = aes(x = p.x, y = p.y, colour = .Tissue), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0(p$axislabel, " 1"), y = paste0(p$axislabel, 
    " 2")) + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("MDS analysis")

ggplot(data = p2) + geom_point(mapping = aes(x = p.x, y = p.y, colour = .Batch), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0(p$axislabel, " 1"), y = paste0(p$axislabel, 
    " 2")) + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("MDS analysis")

geneMatrix <- cpm(y, prior.count = 1, log = TRUE, normalized.lib.sizes = TRUE)

countVar <- apply(geneMatrix, 1, var)

highVar <- order(countVar, decreasing = TRUE)[1:500]
hmDat <- geneMatrix[highVar, ]

pca <- prcomp(t(hmDat), center = TRUE, scale. = TRUE)

loadings <- data.frame(pca$x, .Tissue = metadata$sampleGroup, .Batch = metadata$batch)

ggplot(data = loadings) + geom_point(mapping = aes(x = PC1, y = PC2, colour = .Tissue), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0("PC1, ", round(summary(pca)$importance[2, 
    1] * 100, 2), "% variance"), y = paste0("PC2, ", round(summary(pca)$importance[2, 
    2] * 100, 2), "% variance")) + theme(plot.title = element_text(hjust = 0.5)) + 
    ggtitle("PCA - 500 HVGs")

ggplot(data = loadings) + geom_point(mapping = aes(x = PC1, y = PC2, colour = .Batch), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0("PC1, ", round(summary(pca)$importance[2, 
    1] * 100, 2), "% variance"), y = paste0("PC2, ", round(summary(pca)$importance[2, 
    2] * 100, 2), "% variance")) + theme(plot.title = element_text(hjust = 0.5)) + 
    ggtitle("PCA - 500 HVGs")

tmp <- pheatmap(hmDat, clustering_method = "ward.D2", main = "500 HVGs, all data", 
    show_rownames = FALSE, annotation_col = metadata[, c(1, 3), drop = FALSE], cutree_cols = 4)

dodgy.samples <- grep(paste(c("0412", "1218"), collapse = "|"), colnames(bulk.data), 
    value = TRUE)
metadata$is.dodgy <- ifelse(rownames(metadata) %in% dodgy.samples, "yes", "no")
pheatmap(hmDat, clustering_method = "ward.D2", main = "500 HVGs, all data", show_rownames = FALSE, 
    annotation_col = metadata[, c(1, 3, 4), drop = FALSE], cutree_cols = 4)

Samples originating from patients AHM0412 and AHM1218 are not very good. I will remove them from subsequent analysis.

Double checking after removal of patients

# Names of samples:
samples <- c("0062", "0231", "0412", "0553", "0832", "0848", "0931", "1075", "1218", 
    "1281", "1533", "1566", "1661", "1790")  #RNA-seq data
samples <- samples[!(samples %in% c("0412", "1218"))]  #RNA-seq data
conditions <- c("BE", "GC", "D2", "NE")  # RNA-seq data

files.dir <- "~/Dropbox/Postdoc/2019-12-29_BE2020/RNA-seq/"

# EAC and normal data (tmp, linear scale)
bulk.data <- readRDS(paste0(files.dir, "/EAC-all/count_codingGenes.rds"))
IDs <- grep(paste(samples, collapse = "|"), colnames(bulk.data), value = TRUE)
bulk.data <- bulk.data[, IDs]

# get information about the samples from EAC
metadata <- readRDS(paste0(files.dir, "//EAC-all/samplesAnnotation.rds"))
metadata <- metadata[IDs, 6:8]
group <- factor(metadata$sampleGroup)  #Levels: Barretts Duodenum Gastric Squamous


y <- DGEList(bulk.data, group = group)

# #keep only genes with expression above 1 cpm in at least 3 samples
keep <- rowSums(cpm(y) > 1) >= 3
y <- y[keep, , keep.lib.sizes = FALSE]

# get normalization factors
y <- calcNormFactors(y)


# Start performing diff expression - get sample desing model design <-
# model.matrix(~ batch + group)
design <- model.matrix(~0 + group)
rownames(design) <- colnames(y)

colnames(design) <- levels(group)

# Estimate disparsion
y <- estimateDisp(y, design, robust = TRUE)

plotBCV(y)

# Fit general linear model data
fit <- glmQLFit(y, design, robust = TRUE)

plotQLDisp(fit)

points <- rep(c(0, 1, 2), 2)
colors <- rep(c("blue", "darkgreen", "red"), 2)


p <- plotMDS(y, col = colors[group], pch = points[group], plot = FALSE)


p2 <- data.frame(p$x, p$y, .names = metadata$sampleGroup, .names2 = metadata$batch)


ggplot(data = p2) + geom_point(mapping = aes(x = p.x, y = p.y, colour = .names), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0(p$axislabel, " 1"), y = paste0(p$axislabel, 
    " 2")) + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("MDS analysis")

ggplot(data = p2) + geom_point(mapping = aes(x = p.x, y = p.y, colour = .names2), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0(p$axislabel, " 1"), y = paste0(p$axislabel, 
    " 2")) + theme(plot.title = element_text(hjust = 0.5)) + ggtitle("MDS analysis")

geneMatrix <- cpm(y, prior.count = 1, log = TRUE, normalized.lib.sizes = TRUE)

countVar <- apply(geneMatrix, 1, var)

highVar <- order(countVar, decreasing = TRUE)[1:500]
hmDat <- geneMatrix[highVar, ]

pca <- prcomp(t(hmDat), center = TRUE, scale. = TRUE)

loadings <- data.frame(pca$x, .names = metadata$sampleGroup, .names2 = metadata$batch)

ggplot(data = loadings) + geom_point(mapping = aes(x = PC1, y = PC2, colour = .names), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0("PC1, ", round(summary(pca)$importance[2, 
    1] * 100, 2), "% variance"), y = paste0("PC2, ", round(summary(pca)$importance[2, 
    2] * 100, 2), "% variance")) + theme(plot.title = element_text(hjust = 0.5)) + 
    ggtitle("PCA - 500 HVGs")

ggplot(data = loadings) + geom_point(mapping = aes(x = PC1, y = PC2, colour = .names2), 
    show.legend = TRUE) + theme_bw() + labs(x = paste0("PC1, ", round(summary(pca)$importance[2, 
    1] * 100, 2), "% variance"), y = paste0("PC2, ", round(summary(pca)$importance[2, 
    2] * 100, 2), "% variance")) + theme(plot.title = element_text(hjust = 0.5)) + 
    ggtitle("PCA - 500 HVGs")

tmp <- pheatmap(hmDat, clustering_method = "ward.D2", main = "500 HVGs, all data", 
    show_rownames = FALSE, annotation_col = metadata[, c(1, 3), drop = FALSE], cutree_cols = 4)

dodgy.samples <- grep(paste(c("0412", "1218"), collapse = "|"), colnames(bulk.data), 
    value = TRUE)
metadata$is.dodgy <- ifelse(rownames(metadata) %in% dodgy.samples, "yes", "no")
pheatmap(hmDat, clustering_method = "ward.D2", main = "500 HVGs, all data", show_rownames = FALSE, 
    annotation_col = metadata[, c(1, 3, 4), drop = FALSE], cutree_cols = 4)

pheatmap(hmDat, clustering_method = "complete", main = "500 HVGs, all data, complete clustering", 
    show_rownames = FALSE, annotation_col = metadata[, c(1, 3, 4), drop = FALSE], 
    cutree_cols = 4)

Perform comparison analysis

I increased the treshold for good quality nucleotides. For some reason, the quality of nucletides coming from our analysis if not continous. I only keep nucleotides with phred>34. See the run.txt script in data.path/Counts2 to see how it was calculated.

Read in data

# Names of samples:
samples <- c("0062", "0231", "0412", "0553", "0832", "0848", "0931", "1075", "1218", 
    "1281", "1533", "1566", "1661", "1790")  #RNA-seq data
samples <- samples[!(samples %in% c("0412", "1218"))]  #RNA-seq data
# conditions<-c('BE','GC','D2','NE') # RNA-seq data
conditions <- c("BE", "GC", "NE")  # RNA-seq data

# location of the figures
files.dir <- "~/Dropbox/Postdoc/2019-12-29_BE2020/RNA-seq/"

condition.map <- c(NE = "NE", BE = "BE", GC = "NG")

# Set Conditions
z = 100
vaf = 0.03


# data location
data.path <- "~/Dropbox/Postdoc/2019-12-29_BE2020/Figures/Figure_3/"

# data.path2<-'/home/karolno/Dropbox/Postdoc/2019-06-30_Data/Figure_3E'
# data.path<-'/home/karolno/Dropbox/Postdoc/2019-04-15_Mitotracing/'

# positions to be removed
positions <- c(302:317, 295, 2617, 3107, 13710)



# create list for all the data
all.data <- list()

# Read data all RNA-seq samples all conditions for RNA-seq
for (sample in samples) {
    for (condition in conditions) {
        all.data[[paste0(sample, "_", condition.map[condition], "_RNA")]] <- read.delim(file = paste0(data.path, 
            "/Figure_3E/Counts2/RNA-seq/AHM_", sample, "_", condition, "_Aligned_MT_q.counts"), 
            stringsAsFactors = FALSE)[, c(1:4, 6, 8, 14, 16, 18, 20)]  #keeps information about the count of mutatations
    }
}